home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / studrange.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  198 lines
  1. ; $Id: studrange.pro,v 1.2 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ;  Copyright (c) 1991-1997, Research Systems Inc.  All rights
  4. ;  reserved. Unauthorized reproduction prohibited.
  5.  
  6.  
  7.  function student1_range,Q,V,R,Ifault
  8. ;student_range returns the probability from
  9. ;0 to Q for a studentized range having V
  10. ; degrees of freedom and R samples
  11.  
  12. IFault = 0
  13. if V lt 1.0 or R lt 2. then IFault = 1
  14. if ( Q le 0 or IFault eq 1) THEN $
  15.    return,0 
  16. retval = 0.0
  17.  
  18. QW    = fltarr(30)
  19. VW    = QW 
  20.  
  21. CV    = [.318309886, -.00268132716, .00347222222, .0833333333]
  22. JMin  = 3
  23. JMax  = 15
  24. KMIN  = 7
  25. KMAX =  15
  26. pcutk = .000001
  27. pcutj = .000001
  28. vmax =  120.0
  29. step  = .45
  30. G     = step * R^(-.2)
  31. cvmax = .39894228
  32. GMid  = .5 * ALOG(R)
  33. R1    = R-1
  34. C     = Alog(R*G*cvmax)
  35. if V gt VMAX THEN goto, next20
  36. H     = Step * V ^(-.5)
  37. V2    = .5 * V
  38. CV1   = .193064705
  39. CV2   = .293525326
  40.  
  41. if V EQ 1 THEN C = CV1 ELSE $
  42.  if V eq 2 THEN C = CV2 ELSE  $
  43.     C = sqrt(V2) * CV(0)  $
  44.         / ( 1.0 + (( CV(1) / V2 + CV(2))/V2 +CV(3))/V2)
  45. C    = ALOG(C * R * G * H)
  46.  
  47. next20:
  48.   GStep  = G
  49. QW(0)    = -1.0
  50. QW(JMAX) = -1.0 
  51. PK1      = 1.0
  52. Pk2      = 1.0
  53.  
  54. for k = 1,kmax DO BEGIN
  55.  GStep    = GStep - G
  56. next21:
  57.  Gstep = -  Gstep
  58.  Gk       = GMid + Gstep
  59.  Pk       = 0
  60.  if (pk2 le Pcutk and K gt kmin) THEN goto, next26
  61.  W0 = C - .5 * GK ^2
  62.  PZ = gaussint(GK)
  63.  X =  PZ - gaussint(GK-Q)
  64.  
  65.  if X GT 0 THEN pk = exp( W0 + R1 * ALog(x)) 
  66.  if V GT VMax THEN goto, next26
  67.  jump = -JMax
  68. next22:
  69.  jump  = jump + jmax
  70.  
  71.  for j = 1,jmax Do BEGIN
  72.   JJ = J + jump
  73.   if QW(JJ-1) GT 0 THEN goto,next23
  74.   HJ = H * float(J)
  75.   if ( J LT JMax) THEN QW(JJ) = -1.0
  76.   EHJ = EXP(HJ)
  77.   QW(JJ-1) = Q * EHJ
  78.   VW(JJ-1) = V * (HJ + .5 - .5*EHJ^2)
  79.   
  80. next23:
  81.   PJ = 0
  82.   X = PZ - gaussint(GK-QW(JJ-1))
  83.   If ( X GT 0.0) THEN PJ = EXP(W0 + VW(JJ-1) + R1 * ALOG(x))
  84.   PK = PK + PJ
  85.   if (PJ GT PCUTJ) THEN goto, next24
  86.   if (JJ GT JMIN or K GT KMIN) THEN goto, next25
  87. next24:
  88.  ENDFOR
  89.  
  90. next25:
  91.  H = -H
  92.  if ( H LT 0.0) THEN goto, next22
  93.  
  94. next26:
  95.   retval = retval + PK
  96.  if (K GT KMIN and PK LE PCUTK and Pk1 LE PCutK) THEN $
  97.   goto, done
  98.  PK2 = PK1
  99.  PK1 = PK
  100.  If (Gstep GT 0.0 ) THEN goto,next21
  101. ENDFOR
  102.  
  103. done: return, retval
  104. END     
  105.  
  106.  
  107. function QTRNGO, P, V, R, IFault
  108.  
  109. C1 = .8843
  110. c2 = .2368
  111. c3 = 1.214
  112. c4 = 1.208
  113. c5 = 1.4142
  114. T = gauss(.5 - .5*P)
  115.  
  116. if ( V LT 120.0) THEN T = T + (T^3 + T)/V/4.0
  117. Q = c1 - c2 * T
  118. if ( V LT 120.0) THEN Q = Q - C3/V + C4 * T/V
  119. QTRNGO = T * (Q * ALOG(R-1) + C5)
  120. return,QTRNGO
  121. END
  122.  
  123.  
  124.  
  125.  
  126. function studrange, P ,V ,R
  127. ;+
  128. ; NAME:
  129. ;    STUDRANGE
  130. ;
  131. ; PURPOSE:
  132. ;    Approximate the quantile P for a studentized range distribution 
  133. ;    with V degrees of freedom and R samples. P must be between .9 and .90.
  134. ;
  135. ; CATEGORY:
  136. ;    Statistics.
  137. ;
  138. ; CALLING SEQUENCE:
  139. ;    Result = STUDENT_RANGE(P, V, R)
  140. ;
  141. ; INPUT:
  142. ;    P:    The probability. .9 <= P <= .99.
  143. ;
  144. ;    V:    Degrees of freedom. V >=1.
  145. ;
  146. ;    R:    The number of samples. R > 1.
  147. ;
  148. ; OUTPUT:
  149. ;    The cutoff point of the student_range for probability P, degrees 
  150. ;    of freedom V and sample size R.
  151. ;-
  152.  
  153. NFault = 0
  154. PCut = .000001
  155. JMax = 20 
  156.  
  157. if P  LT .90 or P GT .99 THEN BEGIN
  158.   print,'student_range--- P must be between .9 and .99'
  159.   return, -1
  160. ENDIF
  161.  
  162. if V LT 1 or R LT 2 THEN BEGIN
  163.  print,'student_range, must have V>0 and R >1'
  164.  return, -1
  165. ENDIF
  166.  
  167. Q1 = QTRNGO(P,V,R,NFault)
  168. if Nfault NE 0 THEN goto,done
  169. P1 = student1_range(Q1,V,R,NFault)
  170.  
  171.  if NFault NE 0 THEN goto,done
  172. QTRNG = Q1   
  173. if abs(P1 - P) LT PCut THEN goto,done
  174. if P1 GT P THEN P2 = 1.75 * P -.75 * P1 else $
  175.     P2 = P + .75 * (P - P1) * (1-p)/(1-P1)
  176. if ( P2 LT .8) THEN P2 = .8
  177. if ( P2 GT .995) THEN P2 = .995
  178. Q2 =  QTRNGO(P2,V,R,NFault)
  179.  
  180. if Nfault NE 0 THEN goto,done
  181.  
  182. for j = 2,JMax do BEGIN
  183.  P2 = student1_range(Q2,V,R,NFault)
  184.  if(Nfault ne 0) THEN goto,done
  185.  E1 = P1 - P
  186.  E2 = P2 - P
  187.  QTRNG = (E2 * Q1 - E1 * Q2)/(E2 - E1)
  188.  if abs(E1) LT Abs(E2) THEN goto,next12
  189.  Q1 = Q2
  190.  P1 = P2
  191. next12:
  192.  if abs(P1 - P) LT PCut  THEN goto,done
  193.  Q2 = Qtrng
  194. ENDFOR
  195. done:
  196. return,QTRNG
  197. END
  198.